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Abstract 

In this research, the aim is to develop a repetitive firing stopper mechanism 
using electrical fields exerted on the fiber. The Hodgkin - Huxley nerve 
fiber model is used for modeling the membrane potential behavior. The 
repetitive firing of the nerve fiber can be stopped using approaches based 
on the control theory where the nonlinear Hodgkin - Huxley model is used 
q to achieve this goal. The effects of the electrical field are considered as 

• i-H an additive quantity over the equilibrium potentials of the cell membrane 

channels 
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*sO 1 Introduction 

00 

The repetitive firing of the nerve fibers in the living beings is an important 
subject of the biophysical research since the phenomenon has various meanings. 
It may have a role in several diseases like epilepsy [Pellock et.al (2001)] but 
have also a role in the formation of heart beats where the oscillation frequency 
i-H is in fact represents the heart rate. The Hodgkin - Huxley nerve fiber model 

[Hodgkin et.al (1952)] provides a very useful framework for simulation based 

• neurophysiologic studies where the variation of the membrane potential is mod- 
/\ eled basically as a fourth order nonlinear differential equation. The basic model 

in [Hodgkin et.al (1952)] allows the external current injection as the system in- 
put whereas an additional cell membrane potential can also be considered as an 
input which is the case in [Wang et.al. (2004)]. The electrical fields around the 
nerve cell membranes can induce small electric potentials which may lead to dif- 
ferent behaviors as stated in [Wang et.al. (2004), Kotnik et.al (1998)] however 
the potential induction mechanism can also be used as a repetitive firing control 
mechanism which is the main point of this research. The repetitive firing condi- 
tion on the Hodgkin - Huxley nerve fiber model is detected using a bifurcation 
discovery tool like MATCONT [Dhooge et.al. (2003)] and the washout filter 
[Hassouneh et.al. (2004), Wang et.al. (2000)] and projective control [Medanic 
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(1978)] is integrated to derive a membrane additive voltage profile. This voltage 
profile is converted into an electrical field requirement profile using an approach 
derived from [Kotnik et.al. (1998)]. Using electrical energy in treatment of neu- 
rological disorders is a wide part of neurophysiology research where one of the 
most prominent applications is the deep brain stimulation which is used in the 
treatment of Parkinson's disease [Moro et.al. (2006)] and depression [Mayberg 
et.al. (2005)]. The mechanism is that a device implanted in the brain tissue 
sending electrical pulses to the brain. This can be taught as a current injection 
system whereas in this research the mechanism of electric field induction is used 
to obtain a response from the neurological tissue. The research can be a base 
point for further advancements in electrical field stimulation of the central ner- 
vous system. Using electrical energy in treatment of neurological disorders is a 
wide part of neurophysiology research where one of the most prominent applica- 
tions is the deep brain stimulation which is used in the treatment of Parkinson's 
disease [Moro et.al. (2006)] and depression [Mayberg et.al. (2005)]. The mech- 
anism is that a device implanted in the brain tissue sending electrical pulses 
to the brain. This can be taught as a current injection system whereas in this 
research the mechanism of electric field induction is used to obtain a response 
from the neurological tissue. The research can be a base point for further ad- 
vancements in electrical field stimulation of the central nervous system. In this 
paper, first of all the Hodgkin - Huxley model of nerve fibers together with the 
electrical field effects are reviewed and then the methods of projective control 
theory and washout filters are presented. Next the approach of this research is 
described in detail. 



2 Materials and Methods 



2.1 Hodgkin Huxley Model and Electrical Field Interac- 
tion 

Hodgkin - Huxley model nerve fiber dynamics [Hodgkin et.al. (1952)] is a fourth 
order differential equation derived for modeling the behavior of the membrane 
potential as shown below: 
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The definitions of the parameters are: 

V: the displacement of the cell membrane potential from its resting value in mV 

n : is a dimensional variable which can vary between and 1 representing the 
proportion of activating molecules of the potassium channel 

m : same type of variable as n which represents the proportion of the activating 
molecules of the sodium channel 

h : same type of variable as n which represents the proportion of the inactivating 
molecules of the sodium channel 

I PTi : external current injection in -^4 

cj,i .i cm 

gNa ■ sodium channel conductance in 
gx ■ potassium channel conductance in 

g^: leakage conductance in m ^ due to the chloride and other ions which leads 
to a leaking current flow. Because of that the behaviour is modeled as a con- 
ductance. 

C m : membrane capacitance in < ^ 

Vjvo : sodium channel resting potential in mV 

Vk ■ potassium channel resting potential in mV 

Vl ■ the level of potential where the leakage current reduces to zero 

In the experimentation performed by Alan Hodgkin and his colleagues the nom- 
inal values determined are: 



9Na 




(2.2) 




9l = 0.3 



mS 
cm 2 



V K = -12m V 



V Na = 115m V 



V L = 10.613mF 



In this work the external current injection is set to zero by default. It is to 
be used as an external stimulus for stabilizing the neuron with the washout 
filter. The external effect input to the model in our case is the external voltage 
induction on the cell membrane by the electric fields which is added to (2.1) as: 
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The term Ve is the potential induced on the cell membrane by any external effect 
which is the electrical field in this case. If this is due to a time varying electrical 
field relationship between V E and the exerting electrical field E is obtained from 
the solution of Laplace's partial differential equation which is used by [Kotnik, 
1998] to obtain the following: 

V E = F (s) E (s) Rcos (6) (2.4) 

where E (s) is the electrical field time course defined in the Laplace domain, R 
is the radius of the cell soma and 9 is the directional angle of the electrical field 
exerted upon the nerve fibers. Finally F (s) is the transfer function defined as: 

rp t \ a lS 2 +a 2 s + a 3 

F(s)= M 2 + 6 2S + &3 (2 ' 5) 

with its coefficients: 



ai = 3d\ (Aj (3i? 2 - 3dR + d 2 ) + A m (3dR - d 2 )) 

a 2 = 3d ((A,e Q + X Ei) (3i? 2 - MR + d 2 ) + (X m s + \ e m ) (MR - d 2 )) 

a 3 = 3de Q (e, (3R 2 - 3dR + d 2 ) + e m (3dR - d 2 )) 

h = 2R 3 (A m + 2A ) (A ro + |Ai) +2(R- df (A m - A ) (A, - A m ) 

b 2 = 2R 3 (A, (\e m + e a ) + A m + 2e m + 2e a ) + A ( e< + 2e rn )) +2(R- df 
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The parameters of (2.6) are defined in the following table in reference to [Kotnik, 
1998]: 



Table 1 : The values of electrical parameters for the considered cell membrane 
(values are from Kotnik et.al) 



Parameter 


Definition 


Value 


K 


Cytoplasmic conductivity 


0.3 S-m- L 


Si 


Cytoplasmic permittivity 


7.1 • 10~ iU A - s ■ V~ L ■ m~ L 


A m 


Membrane conductivity 


3- 1CT 7 S -m- 1 




Membrane permittivity 


4A ■ lO' 11 A ■ s ■ V' 1 ■ m' 1 


A 


Extracellular medium conductivity 


3- 1Q~ L S-m' 1 


e 


Extracellular medium permittivity 


7.1 • 10~ iU A - s ■ V~ L ■ m~ L 


R 


Cell soma radius 




d 


Membrane thickness 


5nm 



The relation in (2.4) is derived using admittivity operators. The details of 
derivations are given in the relevant research [Kotnik et.al. (1998)]. In order to 
obtain the level of electric field required for the intended operation the following 
can be derived: 



E ^ = Rlke)* {s)VE (2J) 

where, 

$00= h^+hi+h (2 . 8) 

ais 2 + a 2 s + a 3 

The above can be written since the order of the numerator and denominator in 
(2.5) is same where the inversion does not change the properness of the function. 
For maximum efficiency the directional angle of the field 8 should be zero. 



2.2 Projective Control Theory 

The projective control approach is a linear control methodology that derives an 
output feedback controller from a full state feedback design. The transforma- 
tion from the full state feedback to the output feedback is performed through 
the orthogonal projection which operates on the closed loop eigenspectrum of 
the full state feedback design. Before going into the detail of the projective 
control approach a review of full state feedback control approach is given for 
convenience. Consider a linear system depicted by: 



x = Ax + Bu (2.9) 

where 

xel>el m ,Ae R nxn , B e R nxm 
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If one designs a full state feedback controller as 

u= Kx 

the closed loop dynamics is now 

x = (A — BK) x 

The eigenspectrum of this closed loop dynamics is defined by the equation shown 
below: 



(A-BK)V = VA (2.10) 

whereA is the diagonal matrix with the entries of eigenvalues of the closed loop 
dynamics ( or the eigenvalues of A — BK ) and is the matrix of the corresponding 
eigenvectors. If one is thinking of making a static output feedback from y = Cx 
(C e R qxn , y e R q ) as u = -K D y = -K Q Cx . If this feedback is assumed to 
retain q eigenvalues out of A (denote this as A q and corresponding eigenvectors 
from V as V q ) the following should also be written: 



(A-BKC)V,=V,A, (2.11) 
Combining (2.10) and (2.11) the following can be written: 

(A - BK) V = (A - BK C) V q 
This allows one to solve for 

as: 



K = KV, (CV.r 1 (2.12) 

For the nonlinear systems the model should be linearized using the Jacobian (or 
the Taylor series). 



2.3 Linearization and Bifurcation 



As it is just stated the projective control approach requires linearization of the 
nonlinear model. For the case of a standard nonlinear system representation 
like x = f (x, u) the following equations are written: 



x = Ax + Bu + H(x, u) 
df 



A 
B 



<9x 
du 



(2.13) 
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where x £ R n is the state of the system, u £ R m is the input and H (x, u) is the 
higher order nonlinear terms and Xo is the equilibrium point. The eigenvalues 
of the matrix characterizes the stability of Xo . If the instability of the system 
(2.13) is due to a single pair of complex conjugate eigenvalues the resultant event 
is called as Hopf bifurcation. On the other hand if there is a pair of eigenvalues 
with the same values and opposite signs the resultant issue is called as the saddle 
node bifurcation. For the Hodgkin - Huxley model the repetitive firing in non 
- stimulated (current or other external effects) is a result of the bifurcation 
phenomenon. The bifurcation characteristics of the system can be changed by 
using specially designed controllers. However, many of the controllers lead to a 
change in the position of the equilibrium point. This is not a problem in general 
but if the equilibrium point should not be changed one has to incorporate a 
washout filter which blocks the steady state (DC) inputs to the system. The 
washout filter itself is a linear high - pass filter which can be mathematically 
defined as shown below: 

z = A w z + B w £ 

(2.14) 

tp = A w z + B w £ 

where z £ R p , rp £ R p , £ £ R p , A, w £ R pxp and B w £ R pxp 

If the above dynamics is expressed in the Laplace domain the following result 
can be obtained: 



G (a) = fj*l = A w [sl pxp - A w ] 1 B w + B w (2.15) 

where G (s) is the multi - input and multi - output transfer function from the 
input £ to output ip . It is clearly noted that the following limit tends to zero. 



limG(s)->0 (2.16) 

s— >0 

The above result means that the system in (2.14) blocks the steady state inputs. 
This property is the core part of the control mechanisms aided with a washout 
filter. 



2.4 Hodgkin — Huxley Model and application 

For the Hodgkin - Huxley model in consideration first of all a linearization 



should be done. So the linearization will lead to the linearized model in (2.9 1 
with the matrixA as: 



A = 



dV 
dh 



av 



on on 



On 

o 



dV 
dm 
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1/10 x (1/10 - l/100y)/(exp(l - l/10y) - l) 2 x (1 -n)exp(l - 1/10F) + 1/640 x cxp(-l/80V)n 
^ = -l/10/(exp(5/2 - l/10y) - 1)(1 - m) + 

1/10(5/2 - l/10y)/(exp(5/2 - 1/10V) - 1) 2 (1 - m) x exp(5/2 - I/IW) + 2/9exp(-l/W)m 

dh 

W 
dh 

dn 
dm 

dm 
dh 
dh 

(2.19) 

The terms are given in Section 2.1. The input u and the matrix B is: 



= -7/2000 cxp(-l/20V)(l - h) - l/10/(exp(3 - 1/10V) + ifh x cxp(3 - I/IW) 

- (a n + fi n ) 
= - {a m + Pm) 

= - K + p h ) 



B 



^ (-gN a m 3 h - g K r 






9l) 



u 



Vf, 



(2.20) 



In the next, the repetitive firing condition should be found which is done through 
bifurcation analysis using the MATCONT toolbox of MATLAB ®. For this 
purpose, we will analyze the case where there is a physiological condition on the 
sodium channel. That is the deviation in the equilibrium (rest) potential of the 
channel in which MATCONT provides the following data: 

In order to continue the application, a washout filter should be proposed. The 
input of the washout filter should be a measurable variable which is in this case 
the membrane potential of the nerve fiber. So the washout filter will be of the 
following form: 



z = a w z + b w V 
y = a w z + b w V 



(2.21) 
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Table 2: The repetitive firing condition as a result of deviation in the sodium 
channel 



Parameter & Prop 


V Na = 199.134 


Eigenvalues 


Ai = -5.04011 
A 2 = -0.1257 
A 3 = j0.3958 
A 4 = -j0.3958 


Equilibrium Points 


V = 0.93404938 
n = 0.33208269 
to = 0.059059 
h = 0.5631245 


Type of Bifurcation 


Hopf 



In order to have a stable filter the condition 

a w < 

should be satisfied. For a practical case one can choose a w — —0.01 and b w = 1 
. The resultant filter is augmented into the linearized nerve model in (2.17) - 
(2.21) to obtain: 



A / x / + B f V E 

dV 



dV 

(7m 
dV 
dh 
dV 



dn 







dV 
dm 



chh 
dm 

o 





dV 
dh 




dh 
dh 





B 





V E 



V 
n 
m 
h 

z 



(2.22) 



Numerical representation of the above matrix for the values given in the Table 
2 is now: 



V ' 




- 0.8261 


- 74.9539 


154.0073 


5.3840 







" V ' 




- 0.8261 " 


h 




0.0029 


- 0.1850 













n 







rh 




0.0278 





- 4.0361 










TO 


+ 





h 




- 0.0042 








- 0.1186 







h 







z 




1 











-0.01 




z 








V E 



(2.23) 



The output feedback will be taken from the output of the washout filter so the 
output feedback matrix is C = [ 1 —0.01 ] . In order to proceed, 
a full state feedback controller is necessary for the linear system presented in 
(2.23). This can be obtained using the linear quadratic theory (LQR) and 
the MATLAB's lqr(A, B, Q, R) command directly results a full state feedback 
coefficient matrix K/ where Ve = — K/x . The obtained gain matrix Ky 
minimizes the infinite horizon quadratic gain index shown below: 
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oo 

J = j (x T Qx + u T Ru)dt 



(2.24) 



Taking Q = gl4 X 4 and R = 1 considering the problem of this research the above 
index is rewritten as: 



J = J (q [V 2 +n 2 + m 2 + h 2 ]+V 2 )dt 



(2.25) 



Substituting q = 100 and invoking the necessary command the full state feed- 
back gain coefficient yields: 



K f = [- 10.5426 89.5670 - 133.4628 - 6.4330 - 9.8885] 
The closed loop spectrum (Af — B/K/) for the new case is: 



(2.26) 



5.8641 









3.7002 









1.0170 








0.1849 

- 0.1186 



And the corresponding eigenvectors are: 



0.9937 - 0.9621 0.7096 0.0253 0.0011 

- 0.0003 0.0008 - 0.0024 0.9891 0.0000 

- 0.0057 - 0.0796 0.0065 0.0002 0.0000 
0.0005 - 0.0011 0.0033 0.0016 - 1.0000 

- 0.1122 0.2607 - 0.7046 - 0.1448 - 0.0097 



(2.27) 



(2.28) 



The projective control law will make feedback only from the washout filter thus 
the dimension of the feedback is one and the number of guaranteed retainable 
eigenvalues is also equal to one. So it is best to retain the eigenvalue farthest 
from the imaginary axis which is the one at -8.8641 so the matrix 



is only a single column vector and found as: 



0.99367 

- 0.00032859 

- 0.0057199 
0.00048023 

- 0.11223 



(2.29) 



10 



A)w\MM) 




O SO 100 ISO 200 250 3O0 350 4O0 

Figure 3.1: The repetitive firing response of the bifurcating nerve fiber 

Applying the projection equation in (2.12) gives the feedback gain from the 
washout filter output: 



K Q = KV, (CV^r 1 = - 8.6805 



(2.30) 



And applying the feedback as Af — BfK C yields the output feedback closed 
loop eigenvalues: 



8.8641 

3.1436 

0.0013609 

0.21673 

0.12075 



(2.31) 



So the resultant closed loop is a stable system with the expected eigenvalue at 
-8.8641 in the right place. The control law together with the washout filter is: 



z = -o.oiz + y 

V E = -K (-O.Olz + V) 



(2.32) 



3 Results and Discussion 

The uncontrolled nerve fiber model in (2.1) with the parameters in (2.2) except 
the deviated parameter which is given in Table 2 will yield a repetitive firing 
action potential as shown in Figure 1: 

When the control law of (2.32) is applied through the external voltage input of 
the model in (2.3) the response of the nerve fiber becomes that of the Figure 2. 
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Figure 3.2: The stabilized response of the nerve fiber 

The required level of electrical field can be obtained using the inversion system 
presented in (2.7) with the parameters given in Table 1. Required level of electric 
field is presented in Figure 3. 

As it is expected the level of electric field goes to zero since it stabilized the 
equilibrium condition on the nerve fiber. The steady state portion of the electric 
field is not allowed to pass from the washout filter and thus it stays at zero 
forever after the stabilization. This makes the electric field stimulation of the 
nerve fiber very short. This is a good advantage because of the fact that the level 
of electric field required is quite high as seen from Figure 3. In this research, 
we have presented a stabilization scheme for stopping the repetitive firing of 
the nerve fiber. In order to achieve the goal, the Hodgkin - Huxley model 
of the squid giant axon is taken into consideration. A condition on the sodium 
channel is found where the fiber starts to produce repetitive firing and membrane 
potential is oscillating. After that, the model is linearized and augmented with 
a washout filter to obtain a high pass filtered structure. Finally the projective 
control approach is utilized to produce a static output negative feedback on 
the washout filter output to finish the application. The simulation results show 
that the control laws successfully stabilized the bifurcating nerve fibers with 
a short lasting electric field exertion on the fiber. The disadvantage of the 
method itself is that it requires the bifurcating conditions to be known. If a 
parameter identification method is proposed for the Hodgkin - Huxley model 
by the measurement of the membrane potential for various types of nerve fibers. 
So if the model is properly identified, the research can be beginning point in the 
treatment of diseases based on repetitive firing. In the reverse direction, a stable 
fiber can be forced to repetitive firing condition by using approaches like this. 
For example, if the dendrite and soma membrane potentials can be separately 
measured a pair of pure complex conjugate eigenvalues can be placed so that the 
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Figure 3.3: The required level of electric field for stabilizing the nerve fiber (only 
the first 100 seconds are shown where it stays at zero level after all) 



system is forced into a Hopf bifurcation condition. This will be required since 
in both Hopf and saddle node bifurcations a pair of eigenvalue at the specified 
properties are required. 
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